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Abstract 

The Early Anthropogenic Hypothesis considers that greenhouse gas concentrations should 
have declined during the Holocene in absence of humankind activity, leading to glacial incep- 
tion around the present. It partly relies on the fact that present levels of northern summer 
incoming solar radiation arc close to those that, in the past, preceded a glacial inception 
phenomenon, associated to declines in greenhouse gas concentrations. However, experiments 
with various numerical models of glacial cycles show that next glacial inception may still 
be delayed by several ten thousands of years, even with the assumption of greenhouse gas 
concentration declines during the Holocene. Furthermore, as we show here, conceptual mod- 
els designed to capture the gross dynamics of the climate system as a whole suggest also 
that small disturbances may sometimes cause substantial delays in glacial events, causing 
a fair level of unpredictability on ice age dynamics. This suggests the need of a validated 
mathematical description of the climate system dynamics that allows us to quantify uncer- 
tainties on predictions. Here, it is proposed to organise our knowledge about the physics and 
dynamics of glacial cycles through a Bayesian inference network. Constraints on the physics 
and dynamics of climate can be encapsulated into a stochastic dynamical system. These 
constraints include, in particular, estimates of the sensitivity of the components of climate 
to external forcings, inferred from plans of experiments with large simulators of the atmo- 
sphere, oceans and ice sheets. On the other hand, palaeoclimate observations are accounted 
for through a process of parameter calibration. We discuss promises and challenges raised 
by this programme. 



1 Introduction 



The Early Anthropogenic Hypothesis considers that the natural evolution of climate consisted in 
a decline in greenhouse gas concentrations throughout the Holocene, leading today to conditions 



favourable to accumulation of ice in the Northern Hemisphere (Ruddiman et al. (2011) and 
references therein). This hypothesis supposes an important premise: it is possible to predict 
the slow evolution of climate several millennia ahead. Indeed, suppose that climatologist lived 

8.000 years ago. What the Early Anthropogenic Hypothesis says is that a forecast for the 8,000 
years to come made by this early climatologists would have been a decline in greenhouse gas 
concentrations and ultimately, glacial inception. 

Throughout this paper we will consider that this premise should not be taken for granted. 
The general problem of predicting the evolution of the climate system several millennia ahead 
may be tackled from different perspectives. One method of prediction consists in seeking episodes 
in the past during which the climatic conditions and forcings were analogous to those of the 
Early Holocene, and observe the subsequent climate evolution. The other method consists in 
using models, generally numerical models, which account for known physical constraints on the 
dynamics of the climate system. The two methods are briefly reviewed and discussed in the 
next section of this article. 

Most investigators generally appreciate that the two methods need, in practice, to be com- 
bined. The question, addressed in section 3, is how observations and models may be combined 
in a way that is transparent, optimal, and avoids the dangers of a circular reasoning. The 
methodology proposed here consists in (a) formulating stochastic dynamical systems capturing 
palaeoclimate dynamics. These may be viewed as generators of palaeoclimatic histories, designed 
to encapsulate information inferred from experiments with large numerical models and/or other 
hypotheses about climate system dynamics; (b) calibrate the parameters and initial conditions of 
these dynamical systems on palaeoclimate observations. We explain why the Bayesian statistical 
framework is adapted to this programme. 

2 Traditional approaches to predict the slow evolution of cli- 
mate 

2.1 The analogues method and its limitations 



Among the solutions envisioned by Edwar Lorenz (1969) to forecast weather, was an empiri- 
cal approach also known the 'analogues method'. Considering that similar causes should lead 
to similar effects at least within a finite time horizon, weather predictions may be attempted 
by seeking analogous synoptic situations to today's in the archives of meteorological services. 
Likewise, predicting the natural climate evolution during the Holocene may be addressed by 
identifying in palaeoclimate archives situations for which the state of climate and the forcing 
conditions were similar to those experienced during the Holocene. 

At these times scales climate is influenced by the slow changes in incoming solar radiation, 
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and the latter is a function of the Earth orbital elements (characterised by the eccentricity e 
and the true solar longitude of the perihelion w) and the tilt — or obliquity — of the equator on 
the ecliptic, denoted e. These three elements have all specific and well-known effects on the 
seasonal and spatial distributions of incoming solar radiation (insolation). They also have their 
own periodicities: You never encounter twice a same orbital configuration (see, on this topic, 
the early works of Berger (1977) and Berger (1979])). 

Palaeoclimatologists therefore tend to consider specific measures of insolation, supposed to 
summarise in one quantity the effect of insolation on climate, and look at past times where that 
specific measure of insolation was the same as during the Holocene. However, the choice of an 
'analog' depends on how these measures are constructed (Loutre et al. 2008). 

To illustrate this point consider three measures of insolation classically referred to in the 
tradition of Milankovitch's theory (Figure [TJ: daily mean insolation at the summer solstice 
(21th June), calculated according to Berger (1978a); (b) insolation integrated over the caloric 
summer (the 180 days of highest insolation, assuming a 360-day year) ( jBerger 1978b) — this 
measure is favoured by Ruddiman ( |2007 ) — (c) insolation integrated over all days for which 
insolation exceeds 300 W/m 2 (after a proposal by Huybers and Tziperman (2008)). All are 
calculated here for a latitude of 65° N. It is easily shown that these measures of insolation are 
very well approached by a linear combination of climatic precession (esin(tu)) and obliquity (e); 
they are here ordered by increasing influence of obliquity vs climatic precession. 

Circles on Figure [l] indicate times at which the different insolation measures equal their 
present-day values. Call these 'insolation matches'. Those occurring during a climate environ- 
ment analogous to the present Holocene (arbitrarily defined as a CO2 concentration larger than 
250 ppm and benthic <5 ls O less than 3.8 per mil) are highlighted by vertical bars. 

The timing of insolation matches differ by about 1 ka (thousands of years) depending on 
which insolation is considered: for example the time of insolation match at marine isotopic stage 
(MIS) 11 ranges between —402.9 ka (daily mean at solstice) and —404.5 ka (Huybers' thresholded 
sum). In general, though, insolation levels similar to the present day slightly precede or coincide 
with times during which ice volume markedly increased. Bearing in mind the caveats introduced 
by uncertainties on the chronology of the palaeoclimate records, this observation suggests that 
today's insolation levels are close to those that prevailed at the end of previous interglacial 
periods. However, observe also that present-day insolations are at a minimum on their secular 
course (daily mean solstice) or close to it (season-integrated insolations). This remark invites us 
to be prudent about the prediction of a natural end to the Holocene around the present, again 
assuming no anthropogenic perturbation. 

More crucial differences among insolation measures lie in their dynamics. Climate does not 
respond instantaneously to the astronomical forcing. Ice build-up and ocean alkalinity adjust- 
ments involve time scales of the order of 10,000 years at least. Predicting the evolution of 
climate therefore requires one to consider the history of astronomical forcing over long enough 
a time-window, corresponding to the time needed by the different climate system components 



to adjust dynamically to that specific astronomical history. Ruddiman (2007) and Loutre et al. 
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(2008) noted that the astronomical history preceding —325 ka (marine isotopic stage 9) ad- 



mittedly looks as a reasonable candidate to the last 20 ka by reference to season-integrated 
insolations (consider the 'shapes' of the caloric summer insolation before —325 ka); though, it is 
not straightforward to demonstrate that this really is the 'best possible analog' quantitatively, 
in a way that is convincingly robust and objective (Loutre, pers. comm. and further tests made 
for this article). 

There are further caveats to a straightforward use of insolation analogues to predict climate 
change. First, which insolation (daily mean solstice or seasonal integration) is the best predictor 
of climate change may depend on the level of glaciation (Rohling et al. |2010 ). Second, the 
analogues method supposes as a premise that one single possible climate history may correspond 
to a given history of astronomical forcing. In the formal language of dynamical system theory 



this property is called general synchronisation of climate on astronomical forcing (see Rulkov 



et al. (1995) for a formal definition). It is shown in the next section that this might not be the 
case. Third, climate cycles evidently changed in shape and variance over the last million years 
due to slow changes in the environment. For example, the mid-Pleistocene revolution, when 
glacial cycles shifted from a period of 40 ka to a period of the order of 100 ka is explained by 



Saltzman and Maasch (1990) as a bifurcation related to a slow, tectonically-driven decline in 



background levels of CO2. Clark et al. (2006) review a number of alternative explanations. The 



mid-Briinhes transition about 400 ka ago, after which interglacial periods became characterised 
by higher greenhouse gas concentrations than before this transition, is still unexplained (see 
the review of Tzedakis et al. ( 2009| )). Consequently, one cannot be assured that the levels of 
insolation leading to glacial inception determined on the basis of an event 400 ka ago still hold 
valid today. 



2.2 Using climate simulators to predict the glacial inception 

The other approach to the problem of predicting an ice age proposed by Edwar Lorenz relies on 
a mathematical model of the Earth's climate. In a visionary paper entitled "climate modelling 



as a mathematical problem" (Lorenz, 1970), he expresses the hope that "in the 21st century", 



large mathematical models accounting for "every feature of the atmosphere and its environment" 
as well as ocean circulation and vegetation could be run and predict ice ages. 

The 'mathematical' models envisioned by Lorenz are known today as 'general circulation 
models' or 'global climate models'. They reproduce many features of the climate system visible at 
the global scale (the ocean circulation, modes of variability such as El-Niho) based on equations 
that are expressed at the level of a grid cell. We prefer to call them 'simulators' in order to 
avoid confusion with other meanings of the word 'model'. 

Accounting for physical constraints is naturally considered to strengthen the reliability of 
a prediction, especially if this prediction involves a situation that is unprecedented. Lorenz's 
proposal is therefore a generally well accepted way of dealing with the problems of climate pre- 
dictability at all time scales. However, it has a number of drawbacks that need to be emphasised 
in the present context. An often-heard objection to Lorenz' proposal is computing time. State- 
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of-the-art simulators resolving the dynamics of the ocean and the atmosphere are designed to 
be run on time scales of the order of a century. One experiment of thousand years with such 
a simulator may take months of computing time when it is possible. Experiments over several 
millennia therefore require changes in the simulator design. Modellers have generally chosen to 
degrade, sometimes very drastically, the resolution of the atmosphere and ocean and simplify 
or suppress many of the parametrisations in order to gain computing time. Cloud cover, for 
example, may simply be considered as constant. In the meantime, they attempted to account 
for the dynamics of other components of the climate system, such as the ice sheets and the 
carbon cycle, which react on the longer time scales. The simulators satisfying these criteria are 



known as Earth Models of Intermediate Complexity, or simply EMICS (Claussen et al. 2002). 
EMIC experiments designed to test the early anthropogenic hypothesis were presented soon 



after the original publication of 


Ruddiman 


(2003) 


((Petoukhov et al. 


2000, 


Calov and Ganopolski 



declining CO2 concentrations during the Holocene from about 265 ppm to about 235 ppm (a 
similar decline was imposed on methane concentrations). With that scenario, CLIMBER does 
not accumulate ice any time during the Holocene, while it correctly reproduces the increase in 
ice area during the Eemian / Weschelian transition (115 ka ago) |Claussen et al. (2005). Similar 
experiments were carried out with the LLN-2D simulator of Gallee et al. (1992). Consistent with 
the experiments with CLIMBER, reasonable no-antropogenic scenarios with CO2 declining to 
about 240 ppm cause almost no accumulation before several tens of thousand years in the 
future. Only very drastic declines in CO2 concentrations, such as to reach a CO2 concentration 
of 225 ppm at present and 206 ppm 4 ka in the future, yield some ice accumulation (about 27 m 
of equivalent eustatic sea- level in 22 ka). 

Results of experiments with EMICS may however be questioned because they misrepresent 
(if at all) features of the climate system that may be important to understand and predict its slow 
evolution, such as monsoons, western boundary currents, modes of variability like ENSO and the 
details of the ocean circulation in the Southern Ocean. For these reasons a number of experiments 
were published, based on simulators with more state-of-the-art representations of ocean and 
atmosphere dynamics. These simulators are not designed to compute the accumulation of ice 
but they are useful to indicate spots where snow is likely to accumulate from one year to the next, 
and/or to analyse the conditions propitious for accumulation. They also provide information 
about the sensitivity of the ocean surface and circulation to changes in the astronomical forcing 
and greenhouse gas concentrations. 

In the present special issue, three articles discuss experiments with the Community Climate 



Model in different configurations: Vavrus et al. (2011) examine the influence of topography 
representation on the snow accumulation process and find that accounting for higher resolution 



topography increases the sensitivity of snow accumulation on the external forcing; Kutzbach 



(2011) examine the influence of decreasing greenhouse gas concentrations on sea-surface tem- 
peratures and attempts to quantify the effects on the carbon cycle; they suggest that the feedback 
of the surface cooling on the carbon cycle is substantial enough to accommodate Ruddiman's 
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suggestion of a natural amplification of the natural perturbation Ruddiman (2007). Vettoretti 



and Peltier (2011) uses the ocean- atmosphere version of the Community Climate Model and 



compares the effects of decreasing CO2 concentrations with those of the orbital forcing on snow 
accumulation and the abyssal circulation in the Atlantic. They come to a somewhat challeng- 
ing conclusion for the Early Anthropogenic Hypothesis, that is, astronomical forcing is a more 
important driver of ice accumulation than CO2. Earlier, Vettoretti and Peltier (2004) presented 
a series of eight experiments with the Canadian Climate Centre Atmospheric Model GCM2 
to identify the conditions propitious to year-to-year snow accumulation on ice nucleation sites 
(see also references to closely related work in that article). They found that obliquity greatly 
influences snow accumulation in their model. 

Here, I would like to point to perhaps a more fundamental challenge to Lorenz's vision. Even 
the most sophisticated simulators are an incomplete (truncated) representation of reality. Some 
of the truncated processes are represented with semi-empirical formulae, called parameterisa- 
tions (see Palmer (2005) for a very accessible overview). A fraction of this uncertainty — called 
the parametric uncertainty — can be quantified by considering the sensitivity of predictions on 



parameters involved in the representation of the truncated processes. Pioneering works of Mur- 



phy et al. (2004) on this subject lead to great research activity. However, it will never be 



possible to guarantee that all the physical and biogeochemical processes are correctly and accu- 
rately represented in a numerical model of the climate system. There is structural uncertainty. 
In particular, remember that CO2 is active not only as a greenhouse gas, but also as a fertiliser 
and an ocean acidifier; changes in CO2 concentrations at glacial-interglacial time scales have 
innumerable consequences on life, and hence on climate. 

Saltzman already noted Saltzman (2001)(sect. 5.1) that expecting from a climate simulator 
to reproduce ice ages without any prior information about the timing of glacial-interglacial cycles 
is illusory, so huge would the accuracy needed on the snow accumulation imbalance and carbon 
dioxide fluxes to capture the rate and timing of glacial cycles. Given that glacial cycles cannot 
be simulated ab initio, past observations have to be somehow incorporated into the simulator. 
Very often, the procedure is quite implicit: The uncertain parameters of the climate simulator 
are varied until the simulated climate (or climate changes) is reasonably realistic. This is called 
tuning. One must be clear about the fact that published predictions of the next glacial inception 
are all obtained with climate simulators 'tuned' on past climate observations. For example, the 



LLN-2D simulator, used in Berger and Loutre (2002), Loutre and Berger (2003) and Loutre 



(2003) was tuned to capture the timing and rate of the latest glacial inception, about 115 ka 
ago (Andre Berger, pers. comm.). Confidence in the simulator was gained by the fact that once 
tuned, it reproduces the last glacial-interglacial cycle (Gallee et al. , 1992 ) reasonably well, as well 
as the dynamics of previous interglacial periods (Loutre and Berger , [2003 ) and, admittedly with 



more difficulties, the succession of marine isotopic stages 12-11-10 (400,000 years ago) (Loutre 



2003). 



The problem is that tuning a simulator is a fairly non-transparent and non-optimal way of 
using past climate information for the prediction. It is non-transparent in the sense that it is not 
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always clear which information has been incorporated and which one was effectively predicted (or 
'hindcasted'); it is non-optimal in the sense that only a fraction of the palaeoclimate information 
has been used in the inference process, and that the simulator can never be perfectly calibrated. 
Furthermore, the 'satisfaction' criteria which preside the tuning process are not explicitly defined 
and the consequences of the deviations between simulator results and observations on predictions 
are not explicitly quantified. 



3 Dynamical systems as models of glacial cycles 



3.1 What do we learn from simple dynamical systems? 

The complexity and multitude of mechanisms involved in glacial-interglacial cycles may leave the 
modeller hopeless. Fortunately, this is a fact of nature that very complex systems may exhibit 
fairly regular and structured dynamics. This phenomenon, called 'emergence', is related to the 
fact there are constrains on the system taken as a whole. The dynamics of the system may 
therefore be inferred and understood by identifying the constraints that are relevant at the scale 
of interest without having to consider all the interactions between the individual components of 
the system. 

This argument justifies the interest of scientists for 'conceptual' models of glacial cycles. 
These models do have educational interest and they also have predictive power if they are cor- 



rectly formulated (see examples in Paillard (2001); for more general background about complex 



systems consider the monographs of Haken (2006) and Nicolis and Nicolis (2007)). 

For illustration consider the following minimal conceptual setting. Suppose that the state 
of the climate system may be summarised by two variables. Call them V and D. Variable 
V represents the total continental ice volume. It accumulates variations in time due to the 
astronomical forcing F(t) plus a drift term associated with the variable D, which is defined 
later. Mathematically this reads, if 5V is a variation in V over a fixed time interval 5t: 

5V 



St 



■(D + F(t)+P) 



(la) 



Parameter (3 is introduced later, r is a characteristic time that controls the rate at which ice 
volume grows or decays, and which is tuned on observations. Variable D represents the state 
of a climatic component that responds with a hysteresis behaviour to changes in ice volume. 
The variations in surface pC02 associated with ocean circulation changes may present such a 



behaviour (e.g. (Gildor and Tziperman, 2001 Paillard and Parrenin, 2004)). This is mathemat 



ically modelled as follows: 



SD 
~5t 



a(9D 3 - 3D-V + 1), 



(lb) 



where a is another parameter that controls the rapidity of the transition between the two 
branches of the hysteresis (to understand the principle of this second equation, consider the 
situation V = 1. There are then two possible values of D that would be in equilibrium with V 
(i.e.: ^? = 0): D = 0.57 or D = —0.57. Whichever is approached depends on the system history, 
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i.e.: there is hysteresis). Equations ([T^,b) constitute the dynamical core to several conceptual 
models of ice ages available in the literature, including models by Saltzman and Maasch (1991); 



Gildor and Tziperman (2001) and Paillard and Parrenin (2004). Mathematicians will recognize a 



particular formulation of the celebrated FitzHugh-Nagumo model (e.g. Izhikevich and FitzHugh 



(2006)) originally developed to study neuronal responses to electrical stimuli (the term '+1' 
added at the end of equation ( lb ) plays no other role than constraining V to be almost always 
positive, which is more physical for an ice volume). In the absence of astronomical forcing 
the behaviour of V is mainly dictated by (5. It may converge to a fixed point close to zero 
(if j3 < 1/3), much above zero (/? > 1/3) or oscillate between the two (if —1/3 < f3 < 1/3). 
Within the oscillation regime j3 controls the asymmetry of the spontaneous oscillations and 
a, the sharpness of the shifts between the glaciation and deglaciation regimes. Here we chose 
j3 = 0.2 to capture the notion of slow glacial build-up and fast deglaciation, and a = 3 for 
rapid shifts in D compared to variations in V. The oscillations are also known as relaxation 
oscillations, first documented by Van der Pol in 1921 (Kanamaru, 2007). 

For forcing the system we use a linear combination of climatic precession and obliquity in 
order to capture the notion that summer insolation controls ice ages (F(t) = 30[esinzu + e'], 
where ' is the departure to the mean). With this scaling the forcing amplitude is approximately 
of the order of 1 and thus of the same order as D. There are already two lessons to be drawn 
from this simple experiment. Simulations with this model are shown on Figure [2] First, observe 
that the simulated future climate history is an ice volume decline up to 20 ka into the future, in 
spite of the fact that the model correctly reproduces the increases in benthic trends at the end 
of marine isotopic stages 11 (400 ka ago), 9 (325 ka ago, although the model spuriously produces 
a subsequent drop in ice volume), 7 (200 ka ago) and 5e (115 ka ago, but as for stage 9 the 
subsequent drop towards stage 5c is too strong). Hence, the argument that present ice volume 
should be growing (assuming no anthropogenic perturbation) because it did so during analogous 
situations in past interglacials is not sufficient. It is necessary to demonstrate, in addition, that 
this simple model is inappropriate. Second, observe that D shifts to negative values at the end of 
the deglaciations and, in this model, this shift preconditions glacial inception. This may justify 
why, rather than considering insolation, it has sometimes been decided to seek an analogue to 
the present by considering the timing of terminations, because they are a proxy for the state of 



system variables that cannot be immediately observed. Namely, (EPICA community members 



2004) chose to align the previous deglaciation on the deglaciation that lead to stage 11. Based 



on the observed trends in the <5D, they estimated that glacial inception should not occur before 
more than ten thousand years from now. This alignment was debated on the ground that it 



violates the insolation alignment (Crucifix and Berger, 2006 Ruddiman, 2007) (further informal 
debate and correspondence followed). However, they may be a more fundamental reason why 
any 'alignment', i.e., any blind application of the analogues method, may yield a wrong or at 
least overconfident prediction. 

To see this slightly modify (lb) to add a random number of standard deviation 0.05 every 
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1 ka (assume a normally distributed) : 

5D 



8t 



a{D 6 /3 - D - V + 1) + 0.05aV5t 



(lb 



The random (stochastic) process parameterises unpredictable ocean- atmosphere disturbances 
such as Dansgaard-Oeschger events (cf. Ditlevsen and Ditlevsen ( 2009| )), or volcanic eruptions, 
which impact on climate dynamics. The theory that allows us to represent non- resolved dynam- 
ical processes with stochastic processes was introduced in climatology by Hasselmann (1976). 



Nicolis ( 2005 ) discuss to what extent such parameterisations may also be used to estimate the 



effect of accumulating model errors, and a full volume entirely devoted to stochastic parameter- 



isations was recently published (Palmer and Williams, 2010). 



Together, equations (la) and (lb ) form a stochastic dynamical system. Two trajectories 
simulated with two different realisations of the random numbers are shown on (Figure [3]) . In 
principle the standard deviation of the noise (0.05) is an uncertain quantity that should be 
inferred from realistic computer experiments and / or calibrated on observations. Given that 
the present purpose is merely illustrative we are satisfied with the fact that the noise added to 
the evolution of V and D looks reasonable compared to typical benthic and ice core records. 
Indeed, it is seen that the stochastic disturbances are here too small to compromise the shape of 
glacial cycles and the action of astronomical forcing as a pace-maker, but they are large enough 
to alter the timing of glacial events and shift the succession of glacial cycles. Observe that the 
climatic history shown in black is similar to the deterministic simulation shown in Figure [2} 
The red one undergoes a shift around —400 ka. The length of the simulated interglacials around 
— 400 ka are differ in the two experiments, in spite of the fact that the preceding deglaciations 
are in phase and of similar amplitude. 

A similar phenomenon was is fact visible in the experiments with the conceptual model 
of Paillard (1998), who showed that two future climate scenarios, one with immediate glacial 
inception and one with a long-delayed glacial inception, with only tiny changes in model param- 
eters. It is possible to explain these shifts with arguments of dynamical system theory (article 
in preparation). They are related to the fact that the non-perturbed dynamical system may 
accommodate different plausible climate histories, characterised by different series of deglacia- 
tion timings. At certain times, random perturbations (or small parameter variations) may cause 
a shift from one possible history to another one. The (unsolved) mathematical problem is to 
identify and characterise the moments at which such phase shifts are the most likely to occur. 



The possibility of such shifts contrast with the earlier conclusions of Tziperman et al. (2006). 
Whether phase shifts effectively occurred in the true climatic history is still to be elucidated. 
However, if this phenomenon is confirmed, it implies that the predictive horizon of glacial cycles 
is intrinsically limited and that the analogues method may lead us to be overconfident about 
what can actually be predicted. 
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3.2 Developing dynamical systems of glacial cycles. 



The above example suggests that the gross dynamics of glacial-interglacial cycles may be char- 



acterised by a simple dynamical system. However, Tziperman et al. (2006) previously warned 



us about a form of ambiguity related with the interpretation of such simple dynamical systems 
of ice ages. Indeed, we have not been specific at all about the meaning of the variable D. In 



the model of Paillard and Parrenin ( 2004 ) , the role of D is played by a variable quantifying the 



overturning of the southern ocean circulation. Saltzman and Maasch (1990) choose to introduce 



a third variable in order to distinguish the ocean circulation from the carbon cycle. In their 
model, the non-linear terms appear in the carbon-dioxide equations. More interpretations of the 
meaning of V and D are possible. Ambiguity per se is not necessarily a bad thing. It is useful 
to explore and understand the system dynamics without having to enquire about the details of 
the mechanisms involved. However, it is desirable to attach a more concrete meaning to climatic 
variables in order to explore mechanisms and provide credible predictions outside the range of 
observations. There are two complementary ways of doing this. The first one is to verify that the 
relationships between the different variables are compatible with the information inferred from 



experiments with numerical simulators. For example, eq. (la) encodes a linear relationship 
between the ice mass balance and insolation. How does it compare with climate simulators? 
The second one is to connect the variables with actual palaeoclimate observations. If V and D 
represent ice volume and, say, the North Atlantic Overturning cell, their simulated variations 
should be consistent with the current interpretations of planctonic and benthic foraminifera 
records. 

In this section we focus on how we include information inferred from numerical simulators 
in simple dynamical systems (Figure [4] supports the expose) . Our ambition is to reach a level 
of understanding of the simulator response that is more general than mere predictions of the 
future ice volume evolution obtained with this simulator. To this end we propose to examine 
the bifurcation structure of the simulator in response to varying input parameters. 

To this end, consider again LLN-2D. This model is an implementation of equations of quasi- 
geostrophic motion in the atmosphere and rheological equations for ice-sheet dynamics assembled 
with the purpose of studying glacial cycles. As we noted above the level of sophistication of 
LLN-2D is largely superseded by modern simulators. Nevertheless it is a valid example because 



LLN-2D was used in the debate over the early anthropogenic (Crucifix et al. , 2005, Loutre et al. 



2007) . 



In this example, we attempt to determine the number and stability of steady-states obtained 
with LLN-2D, as a function of the precession parameter. Mathematicians devised very ingenious 
and sophisticated ways of exploring the dynamical structure of large mathematical systems but 
here we follow a very simple and intuitive method, previously used to analyse simulators of 



the ocean-atmosphere Rahmstorf et al. (2005) and ocean-atmosphere-ice-sheet systems Calov 



and Ganopolski (2005). It is known as the 'hysteresis experiment'. This is achieved as follows: 



consider eccentricity of 0.025, obliquity of 23.5°, a CO2 concentration of 210 ppm, and run the 
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simulator until equilibrium with w = 90° (perihelion on the 21st June). Then, slowly vary w 
along the trigonometrical circle (the experiment takes the model equivalent of 400,000 years) 
and record the total continental ice volume as w varies. The resulting curve (Figure [5]) depicts 
an ensemble of quasi-equilibrium points. The first observation is that there is at least one stable 
equilibrium response for all the forcings in the range explored. Secondly, there are two stable 
states for esm{w) within about —0.02 and +0.01. A similar curve was calculated for a CO2 
concentration of 260 ppm. It is similar in shape but the two-stable-state regime now ranges 
from esin(ro) = —0.24 to and the maximum ice volume is 16 10 6 km 3 ). Similar results were 



previously obtained with the CLIMBER-SICOPOLIS simulator (Calov and Ganopolski 2005). 



The bifurcation structure that emerges from LLN-2D was not apparent until we did the 
simulations. We have therefore learned a new bit of information that we are ready to trust, or 
at least to test against the observations. This information may be encoded as a simple equation, 
which is known to have a similar bifurcation structure as the one found in LLN-2D: 

8V = -#(y,C02,e,e,ro,$)(5t (2) 

where 5V is a variation in ice volume over a time-step 5t, <j)y is a function of ice volume, 
astronomical forcing and CO2 and other parameters (summarised by (the full expression is 
given in Appendix 1). 

Equation [2] is a surrogate of the LLN-2D model: it is a simple, fast process that behaves 
dynamically in a way that is similar to LLN-2D. Here, we calibrated the parameter ^ such that 
the surrogate reproduces approximately the behaviour of LLN-2D over the last 800 ka. The 
calibration did not yield a simple best value of "if, but rather a distribution of plausible values 
of the parameter ^ that reproduce reasonably well the model output. The resulting fit is shown 
on Figure |HJ 

By finding a surrogate for LLN-2D we have progressed on two fronts. First, we have distilled 
the vast network of partial differential equations implemented in LLN-2D into a much simpler 
structure. This structure can now be compared with canonical dynamical systems in order to 



comprehend its important dynamical characteristics (Ditlevsen (2009) precisely studied an equa- 
tion similar to ^ in the context of glacial cycles). Secondly, it is possible to efficiently express 
our beliefs and uncertainties about this simulator. Namely, uncertainty about the position of 
tipping points may be expressed as uncertainty on the coefficients (cf. Appendix 1). 

Unfortunately the LLN-2D is a very rough simulator by today's standards. One should 
therefore attempt to correct or augment eq. ([2]) based on the knowledge gained from simulators 
more specifically designed to study the response of the atmosphere and oceans to changes in forc- 
ing and boundary conditions. This is where well-designed experiments with general circulation 
simulators of the ocean and the atmosphere, such as those reviewed in the previous section, may 
be useful. Such plans of experiments may help us to locate points in the parameter space from 
which snow may accumulate from one year to the next. More generally, they should allow us to 
explore the (possibly non-linear) response structure of the atmosphere-ocean system to changes 
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in forcing and boundary conditions. To this end, statisticians have developed two concepts that 
should be useful: 



The experiment design theory. This theory provides practical guidance to efficiently sample 
the space of possible inputs (forcing and parameters) of a numerical model in order to learn 



as much as possible about the simulator with as few experiments as possible (Santner et al. 



2003). 



The Emulator is an interpolator that predicts the response of the climate simulator for other 
parameters than the experiments actually made. It also provides an uncertainty of the 



error associated with the interpolation process (Kennedy and O'Hagan, 2001 Rougier 



et al., 2009) 



Just as the LLN-2D surrogate, the emulator provides us with a fast, efficient, high-level 
description of the general circulation simulator response. It may constitute the basis for non- 
linear parameterisations of, the northern hemisphere snow accumulation balance response to 
obliquity, precession, CO2 and ice volume. In turn, these parameterisations can be included into 
a dynamical model of ice ages. 

At last, it happens that certain processes supposed to relevant at the glacial time scales are 
difficult to simulate reliably with numerical simulators. For example, calculating the response 
of the ocean circulation near the Antarctic ice shelves to changes in sea-level is particularly 
challenging because it involves connexions between processes at very different spacial scales. 



Yet, Paillard and Parrenin (2004) speculated that these circulation changes may play a role on 



carbon cycle dynamics in glacial cycles, and they support their claim with a number qualitative 



arguments. Likewise, Ruddiman (2006) suggests a direct effect of precession on the carbon 



cycle through tropical dynamics and, possible, southern ocean temperatures. In such cases 
where numerical simulation are difficult or considered to be too unreliable, the framework of the 
climate stochastic dynamical model gives us the freedom to formulate hypotheses about the role 
of such processes in the form of simple equations, and then to test them against observations. 



3.3 A statistical framework to accommodate observations 



Many of the ideas expressed above were previously formulated by Barry Saltzman (2001). The 



present programme for investigation and prediction of ice ages proposed here goes further. It 
follows a Bayesian methodology because uncertainties are expressed by means of probability 
distribution functions and they are updated using observations. 



Haslett et al. (2006) depicts the problem quite efficiently: the system defined by eq. (la 



( lb ) may be interpreted as a generator of a priori plausible climate histories. They are a 
priori plausible in the sense that they are generated from a system of equations and probability 
distributions of parameters that express, transparently, our beliefs about the dynamics of the 
climate system. 
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In practice, the calibration of this 'generator' consists in rating these climate histories against 
palaeoclimate observations. This operation consists in attributing more or less weight to the dif- 
ferent possible model parameters, depending on the likelihood of the climate histories that they 
generate. This weighting process constitutes the mechanism by which the a priori uncertainty 
on model parameters is updated by accounting for observations. The more general principle of 
model selection consists, in the presence of two different models, in deciding which one produces 
the climatic histories that are the most likely given the observations. The concept of Bayesian 
model calibration was first applied to the problem of glacial cycle dynamics by |Hargreaves and 



Annan (2002). 



The likelihood is, in that framework, an important quantity. It is formally defined as the 
probability of observing what was observed, if the simulated palaeoclimate history was correct. 
Estimating the likelihood therefore requires to express carefully the connections between mod- 
elled climate histories (the 'model'), the real world, and observations. These connections may be 
visualised by means of a Bayesian graphical network, an illustrative example of which is given 
in Figure [7J The graphic is greatly inspired from that published by Haslett et al. (2006) but it is 
here adapted to the glacial cycles problem. Arrows indicate causal dependencies. Namely, our 
knowledge of benthic 5 18 sampled at a certain depth depends on the 5 ls O that was actually 
preserved at that depth (after, e.g., bioturbation), which itself depends on the age corresponding 
to that depth, as well as on the 5 ls O trapped by the foraminifera at the time corresponding to 
that age. Further up, the generated climate histories depend on the parameters of the climate 
stochastic model and on the external forcing. 

The climate stochastic dynamical system used to generate 'palaeoclimate histories' is thus 
part of a large framework generating potential 'palaeoclimate records', chronologies and sed- 
imentation histories. Haslett and Parnell (2008) provide an algorithm to generate sediment 
chronologies compatible with time markers (radiocarbon dates or, in our case, tephra layers 
and/or magnetic inversions). 

Figure [7] is a simplified diagram; the full network should also include parameters controlling 
the response of proxies on climate. Nonetheless, it is good enough to clarify some of our ideas 
about palaeoclimate how inferences are drawn. For example, palaeoclimate scientists often take 



advantage of astronomical forcing to constrain the chronology of climate records Imbrie et al. 



(1984). The Bayesian graphical network shown here makes it clear that this operation requires 
a climate model (the white ellipse representing climate histories generated by the climate model 
lies between the forcing and observations). 

Bayesian methodologies have become increasingly popular in climate science because they 



provide transparent ways of expressing our uncertainties, modelling choices (Annan and Har- 



greaves, 2006; Rougier, 2007 Rougier and Sexton, 2007 Sanso et al. 2008) and the distance 



estimated between the model and reality (on this specific point see Goldstein and Rougier 



(2009)). 

Calibrating dynamical systems is however a very difficult problem. In general, the number 
of possible climate histories compatible with a model is so big that fairly complex algorithms 
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are needed to reduce as much as possible the number of experiments required to efficiently 
update prior distributions. Palaeoclimates present additional specific challenges related to dating 
uncertainties and the complexity of the environmental factors affecting palaeoclimate archives. 



Hargreaves and Annan ( 2002 ) proposed a straightforward application of an algorithm well known 
by statisticians called 'Markov-Chain Metropolis Hastings'. The algorithm works well as long as 
the calibration period is relatively short and that climate trajectories all cluster around a same 
common history. Unfortunately this may not be the case if local instabilities occur, as in the 
example shown on Figure [3j Furthermore, the algorithm demands large computing resources. 
Crucifix and Rougier ( 2009[ ) then proposed a solution based on a statistical algorithm called 



'particle filter for parameter and state estimation' (Liu and West, 2001). This is a sequential 
filter: it 'sweeps' observations forward in time from the past until the present in order to generate 
a posteriori distributions of parameters and model states at the end of the run. The filter 
seemed to be a powerful solution to the problem posed by the existence of local instabilities. 
Unfortunately (again), further sensitivity studies performed after publication lead us to tone 
down our enthusiasm. The filter performs poorly on the model selection problem because it 
fails to discriminate models on the basis of their long-term dynamics. For example, some of the 
models selected by the filter no longer exhibit glacial-interglacial cycles (unpublished results; 
available on request to the author). 

Hope resides in more advanced Bayesian methods, which combine the Metropolis-Hastings 



strategy with the particle filter (Andrieu et al. 2010). An alternative solution is to calibrate 



the parameters on the basis of invariant summary statistics Wood (2010) using a method known 



an 'Approximate Bayesian Computation' Sisson et al. (2007); Wilkinson (2008). Such statistics 



allow one to characterise a climatic trajectory in a way that is not sensitive to its initial condi- 
tions, nor to the exact timing of climatic events. For example, one may attempt to calibrate the 
dynamical system based on the average duration of a glacial cycle, or in the 'skewness' or 'asym- 



metry' of these cycles (King, 1996). More complex metrics may be envisaged to capture complex 



features associated with the non-linear dynamics of the system (Marwan et al. , 2009 Kantz and 



Schreiber, 2003, and. ref. therein). The use of summary statistics, however, should still respect 
the inference process outlined on Figure [7j An easy mistake is to apply summary statistics on 
climate model outputs and compare them straightforwardly with palaeoclimate data, without 
regard to the data preservation, sampling and possibly interpolation processes. Good summary 
statistics should allow us to efficiently discriminate between climate models and, at the same 
time, they should be robust to hypotheses about the data preservation and sampling processes. 



(1994 


): 


Mudelsee 


(2001 


): 



therein constitute important references in this respect. 



4 Conclusion 

The Pleistocene record suggests that modern insolation is not much above that needed for 
glacial inception. However, the complexity of the climate response to the complex astronomical 
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forcing implies that further theoretical elaborations are needed to transform this statement into a 
reliable prediction about the slow evolution of climate like: 'greenhouse gas should have declined 
during the Holocene, leading to glacial inception'. In particular, there is no perfect insolation 
analogue to the Holocene in the past and, would there be one, it cannot be guaranteed that 
climate would behave exactly the same way as during that hypothetical analogue. We mentioned 
several reasons for this. One is the possibility that small disturbances may, at strategic times, 
delay glacial events by several thousands of years. Locating such 'high-sensitivity' times is an 
attractive challenge. If the Holocene was such a time, any small disturbance, including the 
pre-industrial anthropogenic activity, may have durably inflected the course of climate. 

The dynamical system approach proposed here offers opportunities to tackle this problem 
in the tradition of hypothetic-deductive scientific investigation, subject to the usual criteria of 
model parsimony, prediction skill and accommodation of available observations. We have seen 
that very simple dynamical systems may already roughly reproduce the benthic curve; on the 
other hand, the useful output of more sophisticated numerical simulators may also be captured 
by fairly simply structured surrogates or emulators. These considerations encourage our quest 
for fairly compact and efficient predictive models of glacial-interglacial cycles. They are our main 
hope to credibly accommodate the fairly low amount of information present in palaeoclimate 
records with the flow of information generated by the complex simulators of the climate system. 
Indeed, reduced dynamical models designed to emulate complex numerical simulators contain a 
much smaller number of adjustable parameters and they are way more computationally efficient. 
They are therefore easier to calibrate on palaeoclimate observations. Furthermore, they make it 
possible to infer very general statements about climate predictability based on the large body 
of literature available on dynamical system theory. 

Nevertheless, one should not be deceived by the apparent success our simple model and 
its close parents. Finding a dynamical system that convincingly accommodates more than one 
palaeoclimate data series is not trivial. A convincing model of glacial cycles should arguably 
accommodate observations about the oxygen isotopic composition of foraminifera, greenhouse 
gases, and carbon isotopic data, both in the oceans and in the atmosphere. It should explain 
spectral peaks, stochastic aspects, the asymmetry of cycles and slow changes in the system 
dynamics like those that seemingly occurred over the last million years. So far, this target "is 
still elusive" (Raymo and Huybers, 2008). 

Bayesian inference constitutes a powerful framework to handle knowledge and uncertainties 
associated with palaeoclimate dynamics. This framework will lead us to formulate probabilistic 
statements about the next glacial inception and, more generally, assist our investigations about 
the climate system mechanics and dynamics. 

Such a programme for statistical-dynamical modelling of (palaeo-) climates is in its begin- 
nings and serious challenges are faced. The more general problem of statistical calibration of 
dynamical systems is already quite difficult. Palaeoclimate data bear further difficulties, both 
at the climate system level (intertwining of dynamical structures at different time scales like 
Dansgaard-Oeschger events and glacial cycles) and the observation level (time uncertainties, 
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sediment disturbances, uncertain relationships between climate state and proxy. . . ) It is there- 
fore my opinion that we can't yet express our predictions about natural trends in ice volume 
and greenhouse gases by probabilistic statements, and hence, that natural archives of previous 
interglacials do not falsify natural explanations for increasing greenhouse gas trends during the 
Holocene. 



Appendix : the LLN-2D surrogate 

The form adopted for the surrogate is: 



5V 



-Wo + Pi(V)(V - /9 a ) + Pu(esmw) + (3 e (e) + /3 C InCjSt 



(3) 



where V is ice volume, 5V the ice volume increment over 5t = 1 ka.; C is the CO2 concentration, 
e is a random number with Gaussian distribution, are regression coefficients. V is further 
forced to be positive by imposing a lower bound of — VSt on 5V. 

Priors on regression coefficients are normal as follows (volumes are in 10 6 km 3 , angles 
in radians, carbon dioxide concentrations in ppm and times in ka): (3q ~ A/"(19,3 2 ), j3% ~ 
AA(1.4,0.4 2 )/10 3 , Pi ~ A/"(27,3 2 ), /3 n ~ AA(60,5 2 ), (3 £ ~ AA(72,5 2 ), fa ~ AA(3, l 2 ) (conven- 
tionally, Af(fJ,,a 2 ) is a normal distribution of centre \i and standard deviation a). Posteriors are 



estimated by running a particle filter for state and parameter estimation (Liu and West , 2001 ) on 



an experiment spanning the interval [—800 ka; 0] forced by astronomical forcing (Berger, 1978a) 
and C0 2 data by EPICA flLuethi et akj |2008| > and Vostok jPetit et al.[ |1999[ ). The likelihood 
function used to run the filter assumes that the errors of the surrogate every 1-kyr time step are 
independent and normally distributed with a standard deviation equal to 1. 

The red curve on Figure [6] is a climate simulation with the LLN-2D model, while the blue 
curve is computed using the surrogate with the central estimates for Grey curves are trajec- 
tories obtained with parameters sampled from the posterior distribution. Note that the timing 
of the deglaciations is here constrained by the history of CO2 imposed in these experiments. 
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Time-series generated with a two- variable model, equations (la) and (lb). The 



system is forced by a synthetic mix of obliquity and precession. The variable V 
integrates the forcing with a correction D. D shifts between two states according 
to a hysteresis behaviour controlled by V. In this conceptual model V may be 
interpreted as the total volume of continental ice, and D is often associated to 
the state of the ocean circulation. V is here superimposed to the benthic stack 



by Lisiecki and Raymo (2005), which is appropriately scaled to highlight the 



similarity between the simulation and the palaeoclimate record. 
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Same model as Figure [2] but with a small stochastic forcing added to D (equation 



lb | . The two time-series shown correspond to two realisations of the stochastic 



forcing. Observe the phase shift aroud — 400 ka. It corresponds to temporary loss 

of synchronisation of the system on its (astronomical) forcing 

Route map for the design of a stochastic generator of climate histories. Sta- 
tistical approximations of large numerical simulators (surrogates or emulators) 
are developed on the basis of appropriate plans of experiments. These mathe- 
matical objects are much smaller and tractable than the original simulators, and 
may be coupled with each other and/or with more conceptual models. Expert 
elicitation is required for assessing the uncertainties associated with numerical 
simulators and formulate of prior probability distributions of parameters. The 
resulting stochastic dynamical system may then be used to generate ensembles of 

palaeoclimate histories compatible with this model 

Simulated volume of ice over the northern hemisphere near steady-state as a func- 
tion of precession using the LLN-2D simulator by Gallee et al. (1992), and given 
fixed eccentricity, obliquity and CO2 concentration. The experiment consists in 
varying the precession parameter slowly all along the trigonometric cycle in 400 ka 

of model time 

Simulation of ice volume with the LLN-2D model (red) forced by astronomical 
forcing and variations in CO2 (Petit et al. 1999 Luethi et al. 2008), compared 



with 50 realisations of the emulator (eq. Q). The blue trajectory is obtained with 
the emulator using the central estimates of regression coefficients. All experiments 

are deterministic 

Simplified Bayesian inference network for selection and calibration of stochastic 
dynamical systems of palaeoclimates (cf. Haslett et al. (2006) for more theoret- 
ical background). Rectangles correspond to observations or certain quantities. 
Ellipses are uncertain quantities. Continuous functions of time are symbolically 
marked with (t); functions of depth by (d) and discrete series are marked with 
[i]. The sequence from left to right corresponds to forcing, climate model, natu- 
ral archiving of the climate information (blue), sampling and measure (yellow). 
Modelling consists in defining the relationships symbolised by the arrows and at- 
tributing prior probability distributions to uncertain quantities at the top of the 
chain (e.g.: model parameters). Model validation and selection are judged based 
on whether actual observations are compatible the model (measured by the likeli- 
hood); statistical calibration consists in propagating the information backwards, 
from the observations to parameters in order to make predictions (green). 'Sum- 
mary statistics' are variables such as spectral power or other integrated measures 
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Figure 1. Three measures of insolation covering the period — 800 ka to +100 ka (departure to 



the present) calculated using the Berger (1978a) solution, aligned with a stack of 57 benthic 



foraminifera 5 1S (Lisiecki and Raymo, 2005) and a reconstruction of atmospheric CO2 concen- 



tration (Luethi et al. 2008; Petit et al. 1999). Circles mark insolations equal to the present. 



Times corresponding to CO2 concentrations above 250 ppm and benthic 5 1S below 3.8 per mil 
are highlighted to subjectively mark interglacial conditions. 



24 



V-D model 



l— CM 



CM — i 



Q O 




C 

o 'o 
o 



-700 



I 

-500 



I 

-300 
time (ka) 



CM 



CO 



-100 



100 



Figure 2. Time-series generated with a two- variable model, equations (la) and (lb). The 



system is forced by a synthetic mix of obliquity and precession. The variable V integrates the 
forcing with a correction D. D shifts between two states according to a hysteresis behaviour 
controlled by V. In this conceptual model V may be interpreted as the total volume of continen- 
tal ice, and D is often associated to the state of the ocean circulation. V is here superimposed 



to the benthic stack by Lisiecki and Raymo (2005), which is appropriately scaled to highlight 



the similarity between the simulation and the palaeoclimate record. 
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V-D model with stochastic perturbation 
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Figure 3. Same model as Figure [2] but with a small stochastic forcing added to D (equation 



lb ' ). The two time-series shown correspond to two realisations of the stochastic forcing. Observe 



the phase shift aroud —400 ka. It corresponds to temporary loss of synchronisation of the system 
on its (astronomical) forcing. 
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Strategy to design the 'climate stochastic generator' 
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Figure 4. Route map for the design of a stochastic generator of climate histories. Statistical 
approximations of large numerical simulators (surrogates or emulators) are developed on the 
basis of appropriate plans of experiments. These mathematical objects are much smaller and 
tractable than the original simulators, and may be coupled with each other and/or with more 
conceptual models. Expert elicitation is required for assessing the uncertainties associated with 
numerical simulators and formulate of prior probability distributions of parameters. The re- 
sulting stochastic dynamical system may then be used to generate ensembles of palaeoclimate 
histories compatible with this model. 
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LLN-2D simulator 
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Figure 5. Simulated volume of ice over the northern hemisphere near steady-state as a function 



of precession using the LLN-2D simulator by Gallee et al. (1992), and given fixed eccentricity, 
obliquity and CO2 concentration. The experiment consists in varying the precession parameter 
slowly all along the trigonometric cycle in 400 ka of model time. 
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LLN-2D : surrogate vs simulator 
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Figure 6. Simulation of ice volume with the LLN-2D model (red) forced by astronomical forcing 
and variations in CO2 ( |Petit et al. 1999 Luethi et al. 2008), compared with 50 realisations 
of the emulator (eq. Q)- The blue trajectory is obtained with the emulator using the central 
estimates of regression coefficients. All experiments are deterministic. 
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Bayesian inference network 
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Figure 7. Simplified Bayesian inference network for selection and calibration of stochastic dy- 



namical systems of palaeoclimates (cf . Haslett et al. ( 2006 ) for more theoretical background) . 
Rectangles correspond to observations or certain quantities. Ellipses are uncertain quantities. 
Continuous functions of time are symbolically marked with (t); functions of depth by (d) and 
discrete series are marked with [i]. The sequence from left to right corresponds to forcing, cli- 
mate model, natural archiving of the climate information (blue), sampling and measure (yellow). 
Modelling consists in defining the relationships symbolised by the arrows and attributing prior 
probability distributions to uncertain quantities at the top of the chain (e.g.: model parameters). 
Model validation and selection are judged based on whether actual observations are compatible 
the model (measured by the likelihood); statistical calibration consists in propagating the infor- 
mation backwards, from the observations to parameters in order to make predictions (green). 
'Summary statistics' are variables such as spectral power or other integrated measures which 
summarise several observations. 
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